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The law of angular momentum conservation is naturally linked to the rotational symmetry of the involved 
system. Here we demonstrate theoretically how to break the rotational symmetry of a uniaxial crystal via the 
electro -optic Pockels effect. By numerical method based on asymptotic expansion, we discover the 3D 
structure of polarization singularities in terms of C lines and L surfaces embedded in the emerging light. We 
visualize the controllable dynamics evolution of polarization singularities when undergoing the Pockels 
effect, which behaves just like the binary fission of a prokaryotic cell, i.e., the splitting of C points and fission 
of L lines are animated in analogy with the cleavage of nucleus and division of cytoplasm. We reveal the 
connection of polarization singularity dynamics with the accompanying generation of orbital angular 
momentum sidebands. It is unexpected that although the total angular momentum of light is not conserved, 
the total topological index of C points is conserved. 

The study of polarization of light has a long history. Nowadays polarization has been of broad applications in 
many areas of science and technology, ranging from physics to biology and chemistry \ Recent interest may 
be traced back to the seminal work by Nye who revealed the generic structure of polarization singularities^"^. 
In the general 3D picture, there are two types of polarization singularities: lines along which the polarization is 
purely circular (C lines) and surfaces on which the polarization is purely linear (L surfaces), where the orientation 
and handedness of polarization ellipse are indefinite, respectively. In the context of singular optics, polarization 
singularities are regularly considered as the vector analog of phase singularities or optical vortices in scalar fields^. 
Circular polarizations are in essence associated with spin angular momentum of photons, while optical vortices 
are often studied with twisted photons carrying quantized orbital angular momentum (OAM)^'^. It has been 
demonstrated that manipulating optical beams with vortex lines in the forms of knots or links holds promise for 
future laser technology and optical trapping schemes^"^°. Beyond the uniform polarization in scalar fields, the 
morphology and topology of polarization singularities in vector fields are much richer and subtler, as predicted by 
Dennis and later verified by Flossmann et a/^^'^^. Recent years have witnessed a rapidly growing interest in these 
amazing structures, which are found to appear in the skylight^^, isotropic microchip laser^^, near field nano- 
optics^^, and inhomogeneous anisotropic plates^^. 

Here we report theoretically a rather fascinating phenomenon that in a uniaxial crystal when undergoing 
Pockels effect, the tunable evolution of polarization singularities of emerging light behaves just like the binary 
fission, such as in a prokaryotic cell division^^. Under the control of an externally applied electric field, the splitting 
of C points and fission of L lines can be depicted vividly in analogy with the cleavage of nucleus and division of 
cytoplasm. Polarization singularities in crystals have been indeed studied extensively^ ^"^^ but apparently seldom 
considering the Pockels effect^^'^^ Previously, we demonstrated the capability of using Pockels effect to manip- 
ulate spin and orbital angular momentum in optically active crystals or electro -optic birefrigent crystals^^"^^. In 
contrast, we here aim to show another phenomenon of tunable polarization singularities by electro -optically 
breaking the rotational symmetry of a uniaxial crystal, which therefore lends itself to a flexible and real- 
time manipulation. Our work can also be connected to those reporting the conservation law of angular 
momentum related to the rotational symmetry, such as in isotropic crystals or uniaxial crystals^^'^^"^\ In contrast, 
here the rotational symmetry around the optic axis is slightly broken by the applied electric field, since the specific 
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Figure 1 | The theoretical scheme, (a) The proposed schematic diagram, 
(b) The deformed ellipse of refractive index of KDP when undergoing the 
Pockels effect. 

second-order susceptibilities x^y^ and Xyjcz are activated to respond 
for the deformation of the refractive index ellipsoid of electro -optic 
crystal. Furthermore, we reveal the connection of our observations 
with the accompanying generation of OAM sidebands. OAM side- 
bands have been found intrinsic to reflection due to Goos-Hanchen 
and Imbert-Fedorov shifts^^. But ours is resulted from transverse 
angular anisotropy induced by Pockels effect. Surprisingly, total 
angular momentum of light is not conserved, whereas total topo- 
logical index of C points is conserved. 

Our scheme is sketched in Figure 1. We consider a z-cut uniaxial 
crystal of potassium dihydrogen phosphate (KDP), which is a typical 
electro-optic material belonging to class 42m. The principle refract- 
ive indices yIq = 1.5074, Ue = 1.4669, and the nonvanishing electro- 
optic coefficients = 752 = 8, yes ^ H (in pm/V), respectively^l 
Assume the initial light beam of wavelength 1 = 1.633 fim is a left- 
handed circularly polarized one propagating along the optic axis, 
namely, E(ri ,z = 0)= exp ( — r]_/ 2s^)e+ , where s = 4.59 /im is the 

beam waist and e+ = [l,i] as well as e_ = [1, — /] forms the 

v 2 v2 
circular bases. The incident and exit interfaces of KDP are both 
coated with transparent electrode in order to apply a longitudinal 
electric field Eq =£0^2- When Eq is switched on, a nonlinear polar- 
ization responsible for the Pockels effect is induced^^: P^'^ = 
28oX^^\co, 0): EEo, where £0 is the permittivity of free space, /^^^(co, 
0) the second-order susceptibility tensor related to the Pockels effect, 
and E the light field of frequency cd. Starting from Maxwell equations 
and considering total electric displacement D = £0 + P^^, we have 
the following equation governing the complex amplitude of a prop- 
agating light, 

V^E - V(V-E) + klsr^ + ^o^^^P^^ = 0, (1) 

where ko = 2nlk, 8r = diag(nl,nl,nl) denotes the relative dielectric 
tensor, and ^0 is the magnetic susceptibility in vacuum. 

In this paper, we follow Ciattoni's angular spectrum representa- 
tion method^^ and our recently developed numerical method based 
on asymptotic expansion^^ to calculate the complex vectorial field 
E(rx, z) in the propagating space. As E(rj_, z) is obviously position 
dependent, it is naturally expected that the emerging polarization 
from KDP is also spatially variant. Besides, the polarization should 
be electrically tunable by £0 via the Pockels effect. In general, the 
geometry of polarization ellipse can be completely described in terms 
of four Stokes parameters'^: SQ=ExE^^+EyE^y, Si=ExEl—EyE*, 
S2 = 2Re(£*£^), and S3 =2Im(£*£^). Physically, the orientation of 
major axis and the ellipticity of polarization ellipse are characterized 

by ^ = - arg (Si + /S2) and £ = S3/S0, respectively. In very cross sec- 
tion, L lines are those on which S3 = 0, C points are defined as the 
intercept of the loci Si = 0 and S2 = 0. Besides, C points can be 
classified into some basic types, such as lemon, monstar and star^'^\ 

Results 

We plot in Figure 2 our numerical solution of the 3D polarization 
singularities embedded in the emerging light. It is found that when 



the KDP crystal is undergoing the Pockels effect, C lines and L sur- 
faces attain different morphologies. Figure 2(a) demonstrate the 
simple case when Eq is absent. It looks like a right circular cone with 
the conical surface being L surface and the axis of the cone being C 
line. As the light propagates, the C line is stretching along the pro- 
pagation direction coinciding with optical axis. While £0 is switched 
on, however, the C line is quickly bifurcated into two ones, appearing 
like a pair of compasses with both arms being left-handed circular 
polarization and deflecting from the optic axis. Besides, the L surface 
is then gradually cleaved into two separate sleeves and each encircles 
one C line. By a comparison of Figure 2(a)-2(d), we find that apply- 
ing a larger Eq accelerates the cleavage of L surface. As £0 increases 
from 5.31 kV/cm to 10.62 kV/cm, the cleavage point is brought 
forward from z = 6000 /urn to 3000 /nm or so. 

Figure 2 shows only the frame of the 3D polarization structures. 
One can image that the volume is filled with many polarization 
ellipses of various shapes and orientations. For a better view, we also 
visualize their 2D fine structures in Figure 3, assuming the length of 
KDP crystal is fixed at 2: = 6000 fim. By tuning Eq, we observe a 
fascinating phenomenon that the dynamic evolution of polarization 
singularities when undergoing the Pockels effect just behaves like the 
binary fission of a prokaryote celT^. In Figure 3(a), the KDP crystal is 
pure uniaxial without disturbance {Eq = 0). So the central C-point is 
simply surrounded by one L circle. Here we use "L cell" to describe 
the region that the L line embraces. As £0 is increasing, the splitting of 
C points and fission of L lines can be animated in analogy with the 
cleavage of nucleus and division of membrane in a cell division. 
Specifically, the C point first replicates (like a single DNA), then 
attaches each copy to a different part of L cell. In Figure 3(b) with 
£0 = 5.31 kV/cm, the L cell begins to elongate along x direction, and 
the original and replicate C points are pulled apart to separate poles. 
Then the middle portion of the L cell begins to sink, and a cross wall is 
well developed and formed at £0 — 6.90 kV/cm in Figure 3(c). When 
£0 — 10.62 kV/cm, the L cell has been completely spit into two 
daughters of identical C point, shown in Figure 3(d). Obviously, 
present manipulation on polarization singularities could be flexible 
and fast, since the electro-optic Pockels effect possesses a responsible 
time less than one nanosecond^^. 

It is crucial for us to reveal the underlying reasons that support the 
above interesting features. By analogy between polarization and 
phase singularities, we attribute this to the accompanying angular 
momentum dynamics. In principle, we can express the transverse 
light field in terms of both circular polarizations and spiral harmon- 
ics, 

+ GO 

E±(ri,z)= [E;^irL,zre++Ej-(rj_,zre_]exp{il<l>), (2) 

/= — GO 

GO 

and Wf ^ ~ I 27rr|£;- {r^,z)\^dr (with f\ being the normalized con- 

0 

stant) can thus be interpreted as the weight of each OAM mode. As 
an echo of Figure 3, we show the numerical results of Wf in Figure 4. 
In Figure 4(a), the initial left-handed circularly polarized light with / 
= 0 is partially converted into right-handed one with / = 2 while 
acquires Ih OAM per photon, and therefore conserving total angular 
momentum. This is just the case for a pure uniaxial crystaP^'^°. When 
undergoing Pockels effect, besides the energy transfer from left- to 
right-handed circular component, we find that the energy for each 
circular component is distributed over several neighboring even 
OAM modes, i.e., the OAM sidebands are generated due to mode 
coupling. For left-handed component, the main coupling is to the / = 
±2 modes, with a small efficiency of 1.7%, 2.7% and 4.7% in 
Figure 4(b), 4(c) and 4(d), respectively; while those to other higher 
modes are even weaker, as interfered from insets of Figure 4. For 
right-handed component, the significant mode coupling occurs 
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Figure 2 | The 3D structures of Clines (blue) and 1 surfaces (red). Different electric fields are applied: (a) Eq = 0, (b) Eq = 5.31 kV/cm, (c) Eq = 6.90 kV/ 
cm and (d) Eq = 10.62 kV/cm. All coordinates are in unit of jum. See also the Supplemental information video 1. 
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Figure 3 | The 2D fine structure of C points (blue dots) and L lines (red lines). Green lines denote the streamUnes of major axis of polarization 
ellipse. Under the control of £o> they behave like the binary fission of a prokaryotic cell, where E^ is the same as those in Figure 2. See also the 
supplementary information video 2. 
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Figure 4 | The accompanying generation of OAM sidebands. Green and red bars denote the sidebands for left- and right-handed components, 

respectively, where Eq is the same as those in Figures 2 and 3. Insets are the pseudo sidebands of ^ Wf^ , which "exaggerate" the weak coupling to higher 
even OAM of / = ±4. 



between / = 0 and I = +2. Besides, it is electrically tunable. As has 
been revealed by Angelsky et aP^, there is a relationship between 
topological characteristics of component vortices and polarization 
singularities, namely, C points locate at the vortices of the opposite 
circular component. In our case, the left-handed C points are coming 
from the vortices of right-handed component. For a pure uniaxial 
crystal, only I = +2 vortex exists such that only a left-handed C 
points emerges. It is just the superposition of / = 0 and 1= +2 modes 
that accounts for the formation of a pair of right-handed component 
vortices, and therefore, the formation of a pair of left-handed C 
points. As Eq is increasing, the intensity ratio of / = +2 to / = 0 
decreases from 30%, 16% to 5%, shown in Figure 4(b) to 4(d). As a 
consequence, two vortices are pushed away, so are the C points, see 
Figure 3(b) to 3(d). 

Discussion 

The above generation of OAM sidebands can be well understood 
from the transverse angular anisotropy induced by the Pockels effect. 
According to the refractive index ellipsoid theory^^, we know that, 
with the application of electric field along z direction, the transverse 
isotropy of = riy = n^cannot hold anymore; instead, nx = no — 
nly^^Eo 1 2 and = + nly^^Eo / 2. Thus we can define the quantity, 
A = nx — ny = nly^^Eo, to characterize the broken degree of rotational 
invariance around the optic axis, which is evidently proportional to 
Eq. As is well known, the conservation law of angular momentum is 
naturally linked with the rotational symmetry. So here we expect that 
applying a larger Eq will give rise to a larger nonconservative amount 
of angular momentum. Generally, the angular momentum per 
photon within emerging light can be expressed as a sum of spin 
and orbital parts, namely, /= (/+ + (/- 
By calculation, we obtain the angular momentum change (after 



subtracting H for initial left-handed light): Aj = Oh, AJ = 
-0.5794H, A/=-0.9015H and A/=-1.6008H for Figure 4(a)- 
4(d), respectively, thereby confirming our prediction. 

But, surprisingly enough, the total topological index of C points is 
preserved. In Figure 3(a), the polarization streamlines make up spiral 
branches. As we make a complete circuit clockwise around C point, 
note that the polarization ellipse rotates clockwise through a com- 
plete revolution. Gonsequently, the topological index is Jc ^ + 1» and 
this corresponds to the double degeneracy of the central C point. 
While in Figure 3(b)-3(d), the signed number of turns that the 
streamlines makes around each G point is +1/2, and the number 
of streamlines that terminate on the C point is 1, so each C point is a 
lemon type^\ Therefore, the total topological index is preserved, 
namely, Jc — 1/2 + 1/2 = 1, despite that Eq is changing. A complex 
Stokes field, a = Si + /S2 = \a\exp{i20) {0 is the orientation of 
polarization ellipse), is usefully defined to study the Stokes vortex^°. 
It follows that the index of C points Iq is just half the charge of the 
Stokes vortex. As can be seen from Figure 4, only and always the 
right-handed vortex of charge 2 dominates in the OAM sidebands 
such that Ic ^ 1 remain preserved in each subfigure. 

In conclusion, we have discovered and visualized the interesting 
dynamics evolution of polarization singularities for a light field 
emerging from a uniaxial crystal undergoing the Pockels effect, 
where the splitting of C points and fission of L lines are animated 
in analogy with the cleavage of nucleus and division of cytoplasm in 
the binary fission. Because of the rotational symmetry breaking, we 
find that the total angular momentum of light is not conserved, but 
unexpectedly, the total topological index of G points is conserved. 
We reveled the connection of these findings with the accompanying 
generation of OAM sidebands, as a result of the OAM mode coupling 
induced by the Pockels effect. Our results may supply another per- 
spective of angular momentum conservation law in the context of 
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Figure 5 | The schematic diagram for the domain of integration, (a) The oscillatory integrand cos[z/(/cJ] is plotted in the direction, (b) Top 
view of integrating boundary in two-dimensional Fourier space, where the concentric circles are the oscillating periods of the integrand. In both (a) and 
(b), the red point S denotes the stationary point of function /(/c_^ ky)y the gray area centered at point S denotes the domain of integration, and the shading 
region with red slant lines denotes one oscillating period where a coarse-grained average is performed. 



rotational symmetry breaking, and provide a flexible and fast manip- 
ulation on the polarization singularities. 

Methods 

According to the refractive index ellipsoid theory^^, with the application of longit- 
udinal electric field EqCz, the refractive index ellipsoid of KDP is deformed as, 
^^/'^o + ^^/'^o + ^^/'^e + ^y^-i^oah = 1, where a, b, and c denote the crystalline axes. 
Considering the symmetry of a and b in the ellipsoid equation, we choose a new 
coordinate system x,y, and z, where x and 7 are related to a and ^ by 45° rotation while 
z is parral to c, as illustrated in Fig. 1(b). Then the equation of index ellipsoid 
becomes'^ x^(l/nl-\-y^^Eo)-\-y^(l/nl -yga^o) + ^V^e = 1- % considering the 
nonlinearity induced by the Pockels effect as a perturbation, namely, 763-^0 < <nQ^, 
we know that the principle indices of the new index ellipsoid are approximately: 

= "0 — ^1763^0 /2, ny = no + nly^^^EQ /2 and = rie, where the electro -optically 
induced biaxial anisotropy of rix^ riyj^ Uz just accounts for the broken rotational 
symmetry, as discussed above. Besides, in light of such an approximation, we are able 
to present an efficient numerical solution to Eq. (1) based on the Ciattoni's theory of 
angular spectrum representation and the method of asymptotic expansion. 

The key idea of Ciattoni's angular spectrum representation method is to express the 
light field E(r) = E(r j^, z) using two-dimensional Fourier transformation^^. 



By performing a Fourier integral toward E(ki , z), we can finally write the optical field 
in real space as. 



E(ri ,z) = (i^ki exp (zki •ri)E(ki ,z), 



(3) 



where = xe^ +y^y and k^ =kxhx + kyCy are vectors in the real and momentum 
space, respectively. Along this line, to solve Eq. (1), we expand the light field E(r^ , z) 
in the coordinate system a, b, c by 2-D Fourier transformation, and similarly the 

nonlinear polarization V^^(r_i,z)= | d^k_i exp (/kiTi)P^'^(ki,z), whose compo- 
nents in momentum space are: Pf*^ = — Son'^y^^^QEb, P^^ = — Son'^y^^EQEa, and 
Pf^ = 0. Substituting these transformations into Eq. (1), using the above approxi- 
mation adopted in the index ellipsoid theory, and making a 45° coordinate rotation 
back to X, y, and z, we finally arrive the coupling equations that the light field 
components satisfy. 



dE 

- ik. + KkyEy + \k\{n^ - n\y^^E^llf - k^^^ = 0, 



dE 

- iky + k,kyE, + [kl{n^ + nly^^E^/lf - kl]Ey = 0, 



klnl-kl^ '' dz ^ ' dz^- 



(4) 



(5) 



(6) 



The above equations form a complete description of light propagation in the 
momentum space. Note that the coefficients of Ey in Eq. (4) and E^ in Eq. (5) are both 
electrically tunable with Eq, which play a key role in our flexible manipulation of 
polarization singularities. In our simulation, we consider the incident light is a left- 
handed circularly polarized Gaussian one, namely, E(ri ,z = 0) = exp { — r\/ 2s^)e+ . 

With the boundary condition of Ei (0) = — exp ( — /c^ /2)e + , the second-order 

partial differential equations (4)-(6) can be solved to obtain the solution for E(ki,z). 



E(ri,z)= d^\i^< 



exp {\/YiZ + zki Ti ) H 



exp(V^2^ + 2kiTi) >,(7) 



in which the coefficients /li, /l2> and Ci ~ Cg are given by Eqs. (S14) to (S20) in the 
supplementary materials. This equation is the starting point of all our calculations. 
However, normally, the double Fourier integral in Eq. (7) is difficult to be evaluated 
analytically due to the complex forms of Ci ~ Cg. It can be worked out only for several 
special inputting fields and in the paraxial limit^l We therefore turn to numerical 
method to calculate it in this paper. 

In respect of numerical solutions, the main difficulty is that the integrals in Eq. (7) 
are highly oscillating, where the traditional quadrature methods are invalid^*^''*\ To 
overcome this difficulty, we have introduced in Ref [36] an accurate and numerically 
cheap method based on asymptotic expansion theory. Specifically, one need to cal- 
culate a typical oscillatory integral in Eq. (7), 7(z) = J|^(/cx,/cy) exp[izf(kx,ky)]dkxdky, 
where the function g can be Ci ~ Cg, the function /can be/i ~ /2 with 

/i 2 = a/ — /^i 2 + -kx-\--ky. This double integral can be well approximated by 

z z 

asymptotic expansion^^'^^ (or equally the stationary phase method, see the supple- 
mentary materials), when the value of z is sufficiently large. The essence of the 
asymptotic expansion theory is that the dominated contribution to the highly oscil- 
latory integral /(z) only comes from the small regions in the vicinity of certain critical 
points^^'^\ The most important critical points are the stationary points of function 



dl^df_ 

dkx dK 



f{kx, ky) inside the domain of integration, which satisfy the condition 

The other kinds of critical points are defined on the boundary curve of integration and 
thus are unimportant here, for the boundary of Fourier integral of Eq. (7) is actually at 
infinity. 

Based on this observation, our idea of approximation to Eq. (7) is as follows: 
Instead of integrating over the whole Fourier k space, we only numerically integrate 
over small (and important) regions near the stationary points of functions /i ,2 and 
throw away the contributions from the other oscillatory regions^*^. The part of the 
integral that is thrown away decays exponentially fast as the oscillation increases^\ A 
consequent question is how to choose the domain of integration. We have verified 
that a square integrating boundary yields the smallest error with respect to the 
analytical solution^*^. We illustrate a schematic domain of integration in Fig. 5, where 
the domain is denoted by the gray area centered at the stationary point S. The 
accuracy of this method can be further improved by a coarse-grained average treat- 
ment. For example, we can choose ten equally spaced values d]^ in one oscillating 
period of integrand, e.g., in the slant-line shading regions between two red vertical 
lines in Fig. 5(a) or two red squares in Fig. 5(b), calculate the integral using each d^ as 
the width of integrating domain, and treat their average as the final result. In this way, 
the highly oscillating part far away from the critical point is canceled, leaving the most 
important contribution from the small area close to the critical point. 

Our calculating procedure is as follows: We first numerically find the stationary 
points of functions /i_2, determinate the domains of integration by counting the 
periods of exp(/zfi 2), perform the numerical integration in small regions using Gauss- 
Kronrod quadrature, and finally do average over one oscillating period to get the final 
result. More detailed discussions about the validity and accuracy of this numerical 
method can be found in the supplementary materials. After computing the optical 
fields involving the Pockels effect, we can then explore the dynamics of polarization 
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singularities. It is noted that our method can be extended to any crystals of arbitrary 
point group with the biaxial anisotropy is induced by the Pockels effect, Kerr effect or 
other nonlinear optical effects. 
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